夜間光データ及びSDID法を用いた実証分析
library(tidyverse)
library(sf)
library(terra)
library(openxlsx)
library(mapview)
library(RColorBrewer)
library(osrm)
library(movecost)
library(fs)
library(stringr)
library(synthdid)
mapview::mapviewOptions(fgb=FALSE)
load("data/china.RData")衛星夜間光データとは
ntl1992<-terra::rast(x="data/9828827/Harmonized_DN_NTL_1992_calDMSP.tif")
ntl2018<-terra::rast(x="data/9828827/Harmonized_DN_NTL_2018_simVIIRS.tif")
ntl1992 <- ntl1992 %>%
terra::crop(y=ext(c(123.1,146.9,24.0,45.7)))
ntl2018 <- ntl2018 %>%
terra::crop(y=ext(c(123.1,146.9,24.0,45.7)))1992年日本
terra::plot(ntl1992)
2018年光り輝く日本 若干赤くなってる？
terra::plot(ntl2018)
Synthetic DID法とは
発想
因果推論アプローチの一つであるPotential outcome frameworkにおいて重要なのは 観察しえない反実仮想\(Y_{itD=0}\)をいかにして導き出すかという問題です。
反実仮想というのは例えば、統一した東ドイツ（に当たる地域の）経済\(Y_{itD=1}\)に対する統一しなかった場合の東ドイツ経済\(Y_{itD=０}\)や拓銀が倒産した後の北海道経済（観察可能）に対する拓銀が生存した場合の北海道経済（観察不可能）などです。
今回紹介するSDID法は、本当に簡単に言うと
拓銀が倒産しなかった世界線の北海道を、他のいろんな県を組み合わせて作ろう！
という発想が元になっています。
説明
SDID法は差分の差分法(DID法)と合成コントロール法(SC)を組み合わせた手法です。
DID推定量は以下の通りです （書き方は様々あるが）\[ (\hat{\tau}^{DID}, \hat{\mu}, \hat{\alpha}, \hat{\beta})={\arg\min}_{\mu,\alpha,\beta,\tau} \sum_{i=1}^N \sum_{t=1}^T(Y_{it}-\mu-\alpha_i-\beta_t-D_{it}\tau)^2 \] 各ユニット・時間における結果変数\(Y_{it}\)は定数\(\mu\)、時間効果\(\beta\)、個人効果\(\alpha\)の線形結合で表されます。今回求めたい因果効果をとらえるのは処置変数\(D_{it}\)の係数\(\hat{\tau}^{DID}\)です。このモデルの最小化を行うことでDID推定量は求められます。
一方、合成コントロール法はSynthetic Control Methods for Comparative Case Studies(Abadie et al. 2010) で提案されました。
アイデアとしては、処置を受ける前まで処置群を統制群（処置を受けない群）の重み付け（加重平均）で近似できるような重み \(\hat\omega\) を求めて、その加重平均で処置後の処置群を予測し、その差分を処置効果とするものです。それを表す式は以下の通りです。
\[ \hat{\tau}_{it}^{synth}=Y_{Nt}-\sum_{i=1}^{N-1}\omega_{i}Y_{it} \]
問題の重み\(\hat\omega_{i}\)は処置前までの期間、つまり\(T-1\)期に注目して以下のL2正則化項（罰則付き）二乗誤差最小化問題で求められます。
\[ \hat{\omega}_{i}=\underset{\omega}{\arg\min}\frac{1}{T-1} \sum_{t=1}^{T-1}(\sum_{i=1}^{N-1} \omega_{i} Y_{i t}-Y_{N t})^{2}+ \frac{1}{2}\zeta\|\omega\|_2 \]
…つまり、どちらも処置群と統制群のバイアス（元の違い）と処置前と処置後のバイアス（放っておいても変化するかも）を統制して差分(処置効果）を取りたい訳です。そのなかで、いかに”並行トレンド仮定”という非常につよい仮定を緩められるかという所が問題になっていますね。
ですが、合成コントロール法ユニットごとの重み\(\hat\omega_{i}\)しか使っていません。もちろん、重要な年と重要じゃない年を峻別するような時間の重みも考えられますよね!
そのような時間に対する重み\(\lambda_{t}\)を合わせてDID推定量の中に組み込んでしまったのが今回使用するSynthDID推定量です。(\(\hat\lambda_{t}\)の推定に関しては\(\hat\omega_{i}\)の最小化問題の対象がユニットから時間に変わっただけです。)
\[ (\hat{\tau}^{SDID}, \hat{\mu}, \hat{\alpha}, \hat{\beta})={\arg\min}_{\mu,\alpha,\beta,\tau} (\sum_{i=1}^N \sum_{t=1}^T(Y_{it}-\mu-\alpha_i-\beta_t-D_{it}\tau)^2 \hat{\omega}^{SDID}_i \hat{\lambda}^{SDID}_t) \]
うれしさ
北海道拓殖銀行破綻と道内経済への効果検証
jpn1<-sf::read_sf(dsn="data/JPN_shp/gadm41_JPN_2.shp")
mapview::mapview(x=jpn1)